Examen_licenciatura_137926

library(readr)
library(ggplot2)
library(stringr)
library(dplyr)
library(knitr)
library(tidyr) 
library(class)
library(rpart)
library(lubridate)
library(randomForest)
library(dummies)
set.seed(137926)
pollution_complete <- read_csv("pollution.csv")
row.has.na <- apply(pollution_complete, 1, function(x){any(is.na(x))})
pollution <- pollution_complete[!row.has.na,]
numero_NA <- dim(pollution_complete)[1] -dim(pollution)[1]

pollution<- pollution[order( pollution$year, pollution$month,pollution$day),]
totales<- dim(pollution)[1]

#queremos saber que variables son categoricas y cuales numericas
classes <- sapply(pollution, function(x) class(x))
categoric_cols <- pollution[,which(classes %in% c("character", "factor"))]
numeric_cols <- pollution[, -which(classes %in% c("character", "factor"))]

#cuantos valores unicos --cardinalidad
uniques_num <- sapply(numeric_cols, function(x) unique(x) %>% length())
uniques_num <- as.data.frame(uniques_num)

uniques_cat <- sapply(categoric_cols, function(x) unique(x) %>% length())
uniques_cat <- as.data.frame(uniques_cat)

#valores unicos
uniques_values_num <- sapply(numeric_cols, function(x) unique(x))
uniques_values_cat <- sapply(categoric_cols, function(x) unique(x))

#proporcion valores unicos - uniqueness
uniqueness_num <-  round(uniques_num/totales * 100, 2)
uniqueness_num <- as.data.frame(uniqueness_num)

uniqueness_cat <-  round(uniques_cat/totales * 100, 2)
uniqueness_cat <- as.data.frame(uniqueness_cat)
                                
#checamos si hay vacios 
nan_num <- sapply(numeric_cols, function(x) sum(is.na(x)))
nan_num <- as.data.frame(nan_num)

nan_cat <- sapply(categoric_cols, function(x) sum(is.na(x)))
nan_cat <- as.data.frame(nan_cat)

#para sacar la moda, redondeamos a 2 digitos
my_mode <- function(x){
  if (class(x) %in% c("character", "factor")) {
    table(x) %>%
      which.max() %>%
      names()
  }
  else {
    table(round(x, 2)) %>%
      which.max() %>%
      names()
  }
}


#moda
modes_num <- sapply(numeric_cols, function(x) my_mode(x))
modes_num <- as.data.frame(modes_num)

modes_cat <- sapply(categoric_cols, function(x) my_mode(x))
modes_cat <- as.data.frame(modes_cat)


#min 
mins <- sapply(numeric_cols, function(x) min(x))
mins <- as.data.frame(mins)

#max
maxs <- sapply(numeric_cols, function(x) max(x))
maxs <- as.data.frame(maxs)

#mean
means <- sapply(numeric_cols, function(x) mean(x))
means <- as.data.frame(means)

#median
medians <- sapply(numeric_cols, function(x) median(x))
medians <- as.data.frame(medians)

#1st quantile
first_qtls <- sapply(numeric_cols, function(x) quantile(x)[2]) #revisa la funcion quantile!
first_qtls <- as.data.frame(first_qtls)

#3rd quantile
third_qtls <- sapply(numeric_cols, function(x) quantile(x)[4]) 
third_qtls <- as.data.frame(third_qtls)

#sd
sds <- sapply(numeric_cols, function(x) sd(x))
sds <- as.data.frame(sds)

###generamos nuestra tabla de data profiling
df_categoric <- cbind(uniques_cat, uniqueness_cat, nan_cat, modes_cat)
#no me gusta que tengan en la columna cat... 
names(df_categoric) <- str_replace_all(names(df_categoric), "_cat", "")
names(df_categoric)[2]<-"uniqueness"

df_numeric <- cbind(uniques_num, uniqueness_num, nan_num, mins, maxs, means,
                    sds, medians, modes_num, first_qtls, third_qtls) 
names(df_numeric) <- str_replace_all(names(df_numeric), "_num", "")
names(df_numeric)[2]<-"uniqueness"

Introducción al problema

Una investigación del Colegio del Medio Ambiente de la Universidad de Nanjing relacionó la contaminación con casi un tercio de todas las muertes que se producen en China, ubicando a la polución en el mismo nivel que fumar tabaco como amenaza para la salud pública. El estudio analiza casi 3 millones de muertes en 74 ciudades chinas durante 2013. Los hallazgos revelan que un 31,8% de todas las muertes registradas podrían estar relacionadas con la contaminación, con las grandes de ciudades de Hebei, la provincia que rodea a Beijing, clasificadas entre las peores.Extraído de CNN.

Debido a los grandes problemas que ha causado la contaminación en China y la gran importancia en la salud que tiene, nosotros buscaremos analizar y predecir los niveles de contaminación en la ciudad de Beijing usando métodos de regresión. A lo largo de este trabajo estaremos trabajando con una base de datos que contiene información relacionada a la contaminación de Beijing del 1 de enero de 2010 al 31 de diciembre del 2014, y buscaremos ver que tan bien podemos predecir los niveles de contaminación del año 2014, a partir de los años 2010, 2011, 2012 y 2013 utilizando el método de regresión.

En la siguiente imagen podemos observar la calidad del aire segun el PM2.5:

Calidad del aire según PM2.5

Nuestra base de datos contiene 43824 registros, sin embargo, haciendo un pequeño análisis, vemos que contiene 2073 registros con información faltante. Como nuestro modelo es de regresión nosotros eliminaremos esos registros, quedándonos al final con 41751 registros con información consistente. Además, nuestra base de datos contiene las siguientes columnas:

colnames(pollution)

##  [1] "No"    "year"  "month" "day"   "hour"  "pm2.5" "DEWP"  "TEMP" 
##  [9] "PRES"  "cbwd"  "Iws"   "Is"    "Ir"

Donde: + No: Número de Registro + year: Año en el que se hizo el registro + month: Mes en el que se hizo el registro + day: Día en el que se hizo el registro + hour: Hora en el que se hizo el registro + pm2.5: Concentración de PM2.5 + DEWP: Punto de rocío + TEMP: Temperatura + ** PRES: ** Presión + cbwd: Dirección del aire combinada + iws: Velocidad del aire acumulada + Is: Horas acumuladas de nieve + Ir: Horas acumuladas de lluvia

Data Profiling

Hagamos un análisis de las variables categóricas:

kable(df_categoric, format.args = list(big.mark=",", scientific=F))

	uniques	uniqueness	nan	modes
cbwd	4	0.01	0	SE

La única variable categórica que tenemos es ‘cbwd’ que nos muestra prácticamente la dirección del aire. Sin embargo no nos proporciona mucha información relevante, más que la dirección del viento más predominante es ‘SouthEast’

A continuación analicemos las variables numéricas:

kable(df_numeric, format.args = list(big.mark=",", scientific=F))

	uniques	uniqueness	mins	maxs	means	sds	medians	modes	first_qtls	third_qtls
No	41,751	100.00	25.00	43,824.00	22,279.3439918	12,657.6715189	22,435.00	25	11,466.50	33,261.50
year	5	0.01	2,010.00	2,014.00	2,012.0428014	1.4152337	2,012.00	2013	2,011.00	2,013.00
month	12	0.03	1.00	12.00	6.5133290	3.4540277	7.00	7	4.00	10.00
day	31	0.07	1.00	31.00	15.6862111	8.7856552	16.00	13	8.00	23.00
hour	24	0.06	0.00	23.00	11.5014251	6.9248968	12.00	18	5.00	18.00
pm2.5	581	1.39	0.00	994.00	98.6167996	92.0525716	72.00	16	29.00	137.00
DEWP	69	0.17	-40.00	28.00	1.7513832	14.4342402	2.00	18	-10.00	15.00
TEMP	62	0.15	-19.00	42.00	12.4026730	12.1755584	14.00	24	2.00	23.00
PRES	56	0.13	991.00	1,046.00	1,016.4417858	10.3008171	1,016.00	1014	1,008.00	1,025.00
Iws	2,721	6.52	0.45	565.49	23.8622105	49.6112423	5.37	0.89	1.79	21.91
Is	28	0.07	0.00	27.00	0.0553520	0.7789304	0.00	0	0.00	0.00
Ir	37	0.09	0.00	36.00	0.1948935	1.4182654	0.00	0	0.00	0.00

Es importante mencionar que la variable ‘No’ no da información importante respecto a lo que buscamos predecir, es por ello que no será contemplada para el análisis de regresión. Al observar la variable año, podemos observar que hay un equilibrio en cuanto al número de registros y el año, es decir que se tiene una cantidad bastante similar de registros de cada año. Cabe mencionar que la variable año, día y mes serán variables bastante importantes en nuestro análisis ya que es muy factible que la fecha del año sea muy influente en cuanto a los niveles de contaminación. Es importante mencionar que las variables que también están muy relacionadas con la contaminación son la de temperatura, presión y el punto de rocío. Al observar que el promedio de PM2.5 es de 98, podemos confirmar que la calidad del aire de Beijing es mala, ya que lo pasable es no ser mayor a un PM2.5 de 50.

EDA

En esta gráfica se puede notar que las observaciones de año vs PM2.5 no tienen demasiada relación ya que se tienen todo tipo de muestras y por lo tanto no hay una tendencia obvia.

En la gráfica anterior podemos notar que existe una mayor relación entre Mes y PM2.5 que en la gráfica anterior. De hecho podemos observar que en los primeros meses del año los niveles de contaminación son mayores, y en cambio, a mitad del año tenemos bajos niveles de contaminación. Esta relación Mes-PM2.5 puede que sea de gran interés en el modelo que estaremos analizando posteriormente.

La variable Día no tiene una relación muy obvia respecto a los niveles de PM2.5 en Beijing. Sin embargo, podemos ver que a mitad del mes tenemos un concentraciones más altas de PM2.5. Esta variable puede que sea de interés en los modelos que analizaremos.

En el caso de la variable hora tampoco hay una relación obvia con respecto a los niveles de PM2.5. Puede que esta variable presente baja relación con la variable objetivo.

La variable DEWP, que representa a la temperatura de roció, tiene una relación más visible con los niveles de PM2.5 que las variables anteriores. Lo anterior lo podemos ver que mientras más altos son los niveles de PM2.5, los valores de la temperatura se van acercando más a un valor especifico.

Se puede ver que la Temperatura del ambiente y los niveles de PM2.5 podrían tener alguna especie de relación inversa. Ya que podemos observar que mientras más baja sea la temperatura comienza a crecer los niveles de PM2.5. Será cuestión de checar posteriormente si realmente esta variable está relacionada con la variable objetivo.

La presión parece tener una baja relación con la variable que estamos estudiando pero se puede afirmar que esta relación es positiva, y que en niveles intermedios de presión, tenemos más concentración de altos niveles de PM2.5. Es importante señalar que también es un poco difícil de ver esta relación debido a que la variable de presión es discreta en nuestra base de datos.

Al observar la variable Isw (velocidad del aire) con la variable PM2.5, vemos que es muy probable que estás variables estén relacionadas, definitivamente no es una relación lineal, pero puede que se tenga un comportamiento logarítmico.

Al ser ‘Is’ una variable discreta, es difícil determinar si hay una relación entre estas dos. Sin embargo, vemos que mientras más elevados sean los niveles de Is, vemos que hay bajos niveles de PM2.5.

Las variables Is e Ir presentan el mismo problema, la relación de cada una de ellas con PM2.5 es muy baja pero se podría incrementar si se descartan las observaciones de 0 en las dos variables y es que estas observaciones puede que sean los motivos por los que haya una baja relación.

Matriz de Correlación

A continuación mostraremos la matriz de correlación entre las distintas variables, de tal forma que podamos observar de mejor manera si existe una relación o no con la variable objetivo.

corr_pol <- pollution[, -10]
corr_pol <- corr_pol[,-1]
cor(corr_pol)

##                year         month           day          hour       pm2.5
## year   1.000000e+00 -0.0023829270 -5.251158e-05  0.0002174023 -0.01472089
## month -2.382927e-03  1.0000000000  7.043045e-03 -0.0006955103 -0.02394540
## day   -5.251158e-05  0.0070430448  1.000000e+00  0.0004264350  0.08272270
## hour   2.174023e-04 -0.0006955103  4.264350e-04  1.0000000000 -0.02305084
## pm2.5 -1.472089e-02 -0.0239454000  8.272270e-02 -0.0230508369  1.00000000
## DEWP   7.237139e-03  0.2346075283  3.350855e-02 -0.0217124393  0.17140519
## TEMP   5.558957e-02  0.1722187042  2.286414e-02  0.1495291587 -0.09054952
## PRES  -1.338813e-02 -0.0664022871 -1.047948e-02 -0.0419251061 -0.04725732
## Iws   -6.808283e-02  0.0145205626 -4.932765e-03  0.0587804104 -0.24772651
## Is    -1.955319e-02 -0.0628820855 -3.745665e-02 -0.0024457967  0.01926374
## Ir    -2.630411e-02  0.0388957962 -1.128362e-04 -0.0087237524 -0.05137660
##               DEWP        TEMP        PRES          Iws           Is
## year   0.007237139  0.05558957 -0.01338813 -0.068082829 -0.019553194
## month  0.234607528  0.17221870 -0.06640229  0.014520563 -0.062882085
## day    0.033508552  0.02286414 -0.01047948 -0.004932765 -0.037456651
## hour  -0.021712439  0.14952916 -0.04192511  0.058780410 -0.002445797
## pm2.5  0.171405187 -0.09054952 -0.04725732 -0.247726505  0.019263743
## DEWP   1.000000000  0.82381785 -0.77771385 -0.293062318 -0.034932301
## TEMP   0.823817853  1.00000000 -0.82689320 -0.149536956 -0.094795459
## PRES  -0.777713852 -0.82689320  1.00000000  0.178746430  0.070549295
## Iws   -0.293062318 -0.14953696  0.17874643  1.000000000  0.022641302
## Is    -0.034932301 -0.09479546  0.07054930  0.022641302  1.000000000
## Ir     0.125333362  0.04953414 -0.08052264 -0.009146197 -0.009765282
##                  Ir
## year  -0.0263041074
## month  0.0388957962
## day   -0.0001128362
## hour  -0.0087237524
## pm2.5 -0.0513765991
## DEWP   0.1253333619
## TEMP   0.0495341386
## PRES  -0.0805226355
## Iws   -0.0091461971
## Is    -0.0097652816
## Ir     1.0000000000

Al observar los coeficientes de correlación de nuestra variable objetivo (PM2.5) con las demás variables, vemos que estos están muy cercanos al valor 0, lo cual indica que probablemente no hay una relación lineal con la variable que queremos predecir. Debido a ello, estaremos buscando utilizar regresiones tanto lineales como no lineales para predecir los niveles de PM2.5.

Datos Atípicos

out_year <- boxplot(pollution$year, plot=F)
out_month <- boxplot(pollution$month, plot=F)
out_day <- boxplot(pollution$day, plot=F)
out_hour <- boxplot(pollution$hour, plot=F)
out_DEWP <- boxplot(pollution$DEWP, plot=F)
out_TEMP <- boxplot(pollution$TEMP, plot=F)
out_PRES <- boxplot(pollution$PRES, plot=F)
out_Iws <- boxplot(pollution$Iws, plot=F)
out_Is <- boxplot(pollution$Is, plot=F)
out_Ir <- boxplot(pollution$Ir, plot=F)

El número de outliers de cada variable son:

year: 0
month: 0
day: 0
hour: 0
DEWP: 0
TEMP: 0
PRES: 0
Iws: 4891
Is: 368
Ir: 1739

Después sacar la cantidad de datos atípicos de cada variable se puede ver que solo las variables Iws, Is e Ir los presentan. Consideramos que los datos atípicos de la variable ‘Iws’ podrían ser eliminados, ya que en la gráfica vimos que hay una relación bastante obvia con la variable que queremos predecir, y al eliminarlos puede que mejoremos nuestros modelos de predicción. Sin embargo, los datos atípicos de Is e Ir consideramos que son importantes para el modelo y que por ello no deben ser eliminadas, ya que las horas acumuladas de lluvia o nieve pueden ser un factor importante en la contaminación.

Desarrollo de los modelos

pollution$No <- NULL
pollution <- cbind(pollution, dummy(pollution$cbwd, sep = "_"))
pollution$cbwd <- NULL
colnames(pollution) <- c("year", "month","day","hour","pm2.5","DEWP","TEMP","PRES","Iws", "Is","Ir","cbwd_cv","cbwd_NE","cbwd_NW","cbwd_SE")

train_set <- pollution %>% filter( year >= 2010 & year <= 2013)
test_set <- filter(pollution, pollution$year == 2014)
validation_rows <- sample(dim(train_set %>% filter(year == 2013))[1], 
                          size=round(dim(train_set)[1]*0.085,0),
                          replace=F) + dim(train_set %>% filter( year <= 2012))[1]
validation_set <- train_set[validation_rows,]
train_set <- train_set[-validation_rows,]
validation_set$year <- NULL
train_set$year <- NULL
test_set$year <- NULL

A lo largo de este trabajo estaremos trabajando con distintos conjuntos para entrenar, validar y probar los métodos que usaremos para la regresión. Por ello, el conjunto de entrenamiento serán los registros que tienen fecha del 2010 al año 2013 y este contiene 30279 registros. Por el otro lado, nuestro conjunto de validación estará compuesto por el 8.5% de registros que nuestro conjunto de entrenamiento, sin embargo, este solamente estará compuesto por registros del año 2013; este contiene 2813 registros. Por último, el set de pruebas consistirá en los registros del año 2014 que justamente tiene un tamaño de 8659. Como estamos separando nuestros conjuntos por el año, asumimos de alguna forma que el año no afecta los niveles de contaminación, por ello, eliminaremos esa columna de nuestro set de datos una vez que los conjuntos están separados por año. Sin embargo, debemos tener siempre presente lo mencionado con anterioridad.

Así mismo, la columna categórica cbwd que representa la dirección del aire combinada, será utilizada como una columna Dummie.

Regresión Lineal

A continuación empezaremos por hacer una regresión lineal en nuestro conjunto de datos:

regresion <- lm(pm2.5 ~.,  data=data.frame(train_set))
summary(regresion)

## 
## Call:
## lm(formula = pm2.5 ~ ., data = data.frame(train_set))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -161.76  -50.14  -15.34   31.99  899.27 
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.095e+03  8.196e+01  25.567  < 2e-16 ***
## month       -1.516e-01  1.374e-01  -1.103    0.270    
## day          6.528e-01  5.054e-02  12.915  < 2e-16 ***
## hour         1.197e+00  6.856e-02  17.460  < 2e-16 ***
## DEWP         4.430e+00  6.628e-02  66.841  < 2e-16 ***
## TEMP        -6.529e+00  8.173e-02 -79.884  < 2e-16 ***
## PRES        -1.899e+00  8.032e-02 -23.640  < 2e-16 ***
## Iws         -1.876e-01  9.761e-03 -19.217  < 2e-16 ***
## Is          -2.933e+00  5.460e-01  -5.373 7.82e-08 ***
## Ir          -6.294e+00  2.980e-01 -21.124  < 2e-16 ***
## cbwd_cv     -8.279e-01  1.278e+00  -0.648    0.517    
## cbwd_NE     -2.654e+01  1.579e+00 -16.812  < 2e-16 ***
## cbwd_NW     -2.835e+01  1.229e+00 -23.075  < 2e-16 ***
## cbwd_SE             NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 77.12 on 30266 degrees of freedom
## Multiple R-squared:  0.279,  Adjusted R-squared:  0.2787 
## F-statistic: 976.1 on 12 and 30266 DF,  p-value: < 2.2e-16

En un inicio, intentamos eliminar las variables relacionados con la dirección del aire (cbwd), sin embargo, obtuvimos un peor rendimiento del modelo, primero por que los residuales se alejaban un poco del cero, además de que la R² disminuía. Al tener como coeficiente estimado NA en la variable cbwd_SE podemos afirmar que esa variable tiene una relación completamente lineal con las otras columnas. Podemos notar que la mayoría de las variables proporcionan buena información al modelo de regresión, pero a pesar de esto, la R² de esta regresión lineal es muy baja, ya que es de aproximadamente .279 y por esta razón se puede afirmar que esta regresión lineal no es un buen modelo de regresión porque las variables no explican gran parte de la varianza de la variable objetivo. Lo anterior es es justo lo que observábamos en la matriz de correlación.

A continuación grafiquemos los residuales del modelo:

fitted_values <- regresion$fitted.values
residuals <- regresion$residuals
residuals_df <- data.frame(fitted_values=fitted_values,
                           residuals=residuals)
ggplot(residuals_df, aes(x=fitted_values,
                         y=residuals)) +
  geom_point() +
  geom_hline(yintercept=0)

Al observar la gráfica anterior fácilmente podemos ver que los residuales no tienen una varianza constante, que es justo un requisito para que una regresión lineal sea buena. De hecho, no siguen una distribución normal como observamos en el siguiente histograma:

ggplot(residuals_df, aes(x=residuals)) +   geom_histogram()

Su distribución es más parecida a una chi-cuadrada pero esta no es la distribución que buscamos en los errores de una regresión lineal. Lo podemos visualizar más fácil con el QQ-Plot:

#produce una grafica con el QQ de una normal
qqnorm(residuals_df$residuals, 
       ylab="Residuals", main="")
#se agrega la linea de 
qqline(residuals_df$residuals)

Con esto se puede confirmar lo que el histograma mostraba, que los residuos no siguen una distribución normal y es por esto que la regresión lineal utilizada no puede considerarse como un buen modelo para lo que intentamos hacer.

A pesar de que ya vimos que no es un buen modelo, tratemos de predecir los valores de PM2.5 del set de validación:

prediction <- predict(regresion, data.frame(validation_set))
residuos <- prediction-validation_set$pm2.5

df_residual <- data.frame(fitted_values=prediction,
                 residuals=residuos)


ggplot(df_residual, aes(x=fitted_values, y=residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color="red") +
  theme_bw()

ggplot(df_residual, aes(x=residuals)) +   geom_histogram()

MSE <- mean(residuos^2)
MRSE <- sqrt(MSE)
y <- validation_set$pm2.5
R2 <- 1 - sum((y-prediction)^2)/sum((y-mean(y))^2)

qqnorm(df_residual$residuals, 
       ylab="Residuals", main="")
#se agrega la linea de 
qqline(df_residual$residuals)

Como vimos en el set de entrenamiento, se ve claramente que los errores no tienen una distribución normal. Así mismo en el set de validación tenemos una R² de 0.2808162 que es muy baja, así como un MSE de 7231.2048527 y un MSRE de 85.0364913.

Transformación de variable Iws

Como se vio al realizar el EDA de los datos, existe una variable que tiene gran relación con la variable que estamos analizando, pero el problema es que no es una relación lineal, estamos hablando de la variable Iws, que representa la dirección del aire. Es por esta razón que se puede considerar como una alternativa el transformar esta variable para que la relación que presente con PM2.5 sea lineal y al momento de utilizarla en una regresión lineal se pueda incrementar el valor de R².

La transformación que vamos a utilizar es sacar el logaritmo natural de Iws y ver si la relación con Pm2.5 ya es lineal:

train_set_transform <- data.frame(train_set)
validation_set_transform <- data.frame(validation_set)
train_set_transform$Iws <- log(train_set_transform$Iws)
validation_set_transform$Iws <- log(validation_set_transform$Iws)
ggplot(train_set_transform,aes(x = pm2.5, y=Iws ))+geom_point()+ggtitle("PM2.5 vs ln(Isw)")

Después de realizar la tansformación de la variable se puede ver que la relación ya tiene mas pinta de que es lineal, por lo cual se realizará una nueva regresión pero ya no se utilizará la variable Iws normal, sino la variable que se obtuvo del logaritmo de Iws.

regresion_transf <- lm(pm2.5 ~.,  data=data.frame(train_set_transform))
summary(regresion_transf)

## 
## Call:
## lm(formula = pm2.5 ~ ., data = data.frame(train_set_transform))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -164.95  -49.30  -14.32   31.48  896.12 
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2058.47929   81.42194  25.282  < 2e-16 ***
## month         -0.31650    0.13629  -2.322   0.0202 *  
## day            0.65986    0.05020  13.144  < 2e-16 ***
## hour           1.26311    0.06817  18.530  < 2e-16 ***
## DEWP           4.16023    0.06713  61.971  < 2e-16 ***
## TEMP          -6.16345    0.08302 -74.241  < 2e-16 ***
## PRES          -1.84667    0.07981 -23.137  < 2e-16 ***
## Iws          -10.36989    0.36906 -28.098  < 2e-16 ***
## Is            -2.32000    0.54309  -4.272 1.94e-05 ***
## Ir            -6.02758    0.29623 -20.348  < 2e-16 ***
## cbwd_cv      -16.29908    1.42038 -11.475  < 2e-16 ***
## cbwd_NE      -34.14681    1.60380 -21.291  < 2e-16 ***
## cbwd_NW      -28.55126    1.20809 -23.633  < 2e-16 ***
## cbwd_SE             NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 76.6 on 30266 degrees of freedom
## Multiple R-squared:  0.2888, Adjusted R-squared:  0.2885 
## F-statistic:  1024 on 12 and 30266 DF,  p-value: < 2.2e-16

Podemos observar que gracias a esta transformación obtenemos un aumento en R², lo cual es bueno. Sin embargo, no es lo suficientemente alto como para decir que es un buen modelo. Lo anterior lo podemos observar en las predicciones en el set de validación:

prediction <- predict(regresion_transf, data.frame(validation_set_transform))
residuos <- prediction-validation_set_transform$pm2.5

df_residual <- data.frame(fitted_values=prediction,
                 residuals=residuos)


ggplot(df_residual, aes(x=fitted_values, y=residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color="red") +
  theme_bw()

ggplot(df_residual, aes(x=residuals)) +   geom_histogram()

MSE <- mean(residuos^2)
MRSE <- sqrt(MSE)
y <- validation_set$pm2.5
R2 <- 1 - sum((y-prediction)^2)/sum((y-mean(y))^2)

qqnorm(df_residual$residuals, 
       ylab="Residuals", main="")
#se agrega la linea de 
qqline(df_residual$residuals)

Así mismo en el set de validación tenemos una R² de 0.2979386 que es muy baja pero más grande que el primer modelo de regresión lineal, así como un MSE de 7059.0438413 y un MSRE de 84.0181161. Después de hacer la regresión se puede ver que, al igual que la primera regresión, los residuos obtenidos no cumplen con los supuesto que requiere una regresión lineal. Es decir, no tienen una distribución normal y al graficarlos parecen tener cierta tendencia. Por estas razones se puede decir que los residuos presentan heterocedasticidad.

Arboles de regresión

No todo conjunto de datos puede ser aproximado por una regresión lineal, ya que el comportamiento del conjunto de datos puede ser no lineal, como vimos en la sección anterior. Debido a ello optaremos por usar regresión utilizando árboles, que como sabemos es una alternativa para modelos no lineales. A continuación realizaremos nuestra regresión a partir de un mínimo de 5 observaciones por hoja; a continuación veamos graficas de R^2 y de los errores relativos de acuerdo al número de splits del algoritmo:

reg_tree <- rpart(train_set$pm2.5 ~., data=data.frame(train_set), 
                  control=list(minsplit=5),
                  method="anova")
summary(reg_tree)

## Call:
## rpart(formula = train_set$pm2.5 ~ ., data = data.frame(train_set), 
##     method = "anova", control = list(minsplit = 5))
##   n= 30279 
## 
##            CP nsplit rel error    xerror       xstd
## 1  0.08948337      0 1.0000000 1.0000372 0.01519335
## 2  0.06767640      1 0.9105166 0.9129470 0.01471503
## 3  0.02792135      2 0.8428402 0.8494411 0.01342783
## 4  0.02122454      3 0.8149189 0.8201684 0.01346545
## 5  0.01977166      5 0.7724698 0.7693631 0.01290558
## 6  0.01763057      6 0.7526981 0.7577532 0.01275528
## 7  0.01469707      7 0.7350676 0.7414791 0.01220543
## 8  0.01118905      8 0.7203705 0.7340502 0.01207282
## 9  0.01018529      9 0.7091814 0.7155016 0.01192353
## 10 0.01000000     10 0.6989961 0.7095932 0.01190004
## 
## Variable importance
##    DEWP    TEMP   month     Iws    PRES cbwd_NW      Is cbwd_cv 
##      34      20      15      12       9       7       1       1 
## 
## Node number 1: 30279 observations,    complexity param=0.08948337
##   mean=98.53189, MSE=8246.557 
##   left son=2 (6981 obs) right son=3 (23298 obs)
##   Primary splits:
##       DEWP    < -11.5  to the left,  improve=0.08948337, (0 missing)
##       Iws     < 13.845 to the right, improve=0.07352004, (0 missing)
##       cbwd_NW < 0.5    to the right, improve=0.05529398, (0 missing)
##       cbwd_cv < 0.5    to the left,  improve=0.02595287, (0 missing)
##       PRES    < 1032.5 to the right, improve=0.02565387, (0 missing)
##   Surrogate splits:
##       PRES  < 1028.5 to the right, agree=0.834, adj=0.278, (0 split)
##       TEMP  < -3.5   to the left,  agree=0.829, adj=0.257, (0 split)
##       month < 2.5    to the left,  agree=0.819, adj=0.215, (0 split)
##       Iws   < 94.555 to the right, agree=0.790, adj=0.089, (0 split)
## 
## Node number 2: 6981 observations,    complexity param=0.01469707
##   mean=48.90603, MSE=3506.546 
##   left son=4 (4064 obs) right son=5 (2917 obs)
##   Primary splits:
##       Iws     < 8.72   to the right, improve=0.14991590, (0 missing)
##       DEWP    < -15.5  to the left,  improve=0.08052178, (0 missing)
##       TEMP    < -5.5   to the right, improve=0.05496506, (0 missing)
##       cbwd_NW < 0.5    to the right, improve=0.05359431, (0 missing)
##       cbwd_cv < 0.5    to the left,  improve=0.04992375, (0 missing)
##   Surrogate splits:
##       cbwd_NW < 0.5    to the right, agree=0.744, adj=0.388, (0 split)
##       cbwd_cv < 0.5    to the left,  agree=0.705, adj=0.294, (0 split)
##       TEMP    < -5.5   to the right, agree=0.615, adj=0.079, (0 split)
##       cbwd_SE < 0.5    to the left,  agree=0.605, adj=0.054, (0 split)
##       cbwd_NE < 0.5    to the left,  agree=0.599, adj=0.040, (0 split)
## 
## Node number 3: 23298 observations,    complexity param=0.0676764
##   mean=113.4018, MSE=8707.809 
##   left son=6 (18764 obs) right son=7 (4534 obs)
##   Primary splits:
##       TEMP    < 4.5    to the right, improve=0.08329596, (0 missing)
##       month   < 2.5    to the right, improve=0.07361915, (0 missing)
##       Iws     < 7.145  to the right, improve=0.04431538, (0 missing)
##       cbwd_NW < 0.5    to the right, improve=0.02104531, (0 missing)
##       PRES    < 1017.5 to the left,  improve=0.01588710, (0 missing)
##   Surrogate splits:
##       month < 3.5    to the right, agree=0.874, adj=0.351, (0 split)
##       DEWP  < -3.5   to the right, agree=0.864, adj=0.302, (0 split)
##       PRES  < 1023.5 to the left,  agree=0.851, adj=0.236, (0 split)
##       Is    < 0.5    to the left,  agree=0.818, adj=0.066, (0 split)
## 
## Node number 4: 4064 observations
##   mean=29.4813, MSE=975.6753 
## 
## Node number 5: 2917 observations
##   mean=75.9688, MSE=5774.505 
## 
## Node number 6: 18764 observations,    complexity param=0.02792135
##   mean=100.1631, MSE=6430.408 
##   left son=12 (4473 obs) right son=13 (14291 obs)
##   Primary splits:
##       cbwd_NW < 0.5    to the right, improve=0.05778121, (0 missing)
##       Iws     < 3.125  to the right, improve=0.03485716, (0 missing)
##       DEWP    < 10.5   to the left,  improve=0.03428548, (0 missing)
##       month   < 3.5    to the right, improve=0.02271331, (0 missing)
##       cbwd_SE < 0.5    to the left,  improve=0.01720398, (0 missing)
##   Surrogate splits:
##       Iws   < 99.215 to the right, agree=0.779, adj=0.075, (0 split)
##       DEWP  < -10.5  to the left,  agree=0.764, adj=0.008, (0 split)
##       PRES  < 1029.5 to the right, agree=0.762, adj=0.002, (0 split)
##       TEMP  < 38.5   to the right, agree=0.762, adj=0.000, (0 split)
##       month < 11.5   to the right, agree=0.762, adj=0.000, (0 split)
## 
## Node number 7: 4534 observations,    complexity param=0.01977166
##   mean=168.1901, MSE=14405.76 
##   left son=14 (697 obs) right son=15 (3837 obs)
##   Primary splits:
##       Iws   < 23.245 to the right, improve=0.07558569, (0 missing)
##       month < 1.5    to the right, improve=0.06669446, (0 missing)
##       DEWP  < -7.5   to the left,  improve=0.03740316, (0 missing)
##       PRES  < 1026.5 to the right, improve=0.02833582, (0 missing)
##       day   < 16.5   to the left,  improve=0.02426763, (0 missing)
##   Surrogate splits:
##       Ir   < 10     to the right, agree=0.851, adj=0.030, (0 split)
##       Is   < 10.5   to the right, agree=0.850, adj=0.022, (0 split)
##       PRES < 1039.5 to the right, agree=0.846, adj=0.001, (0 split)
## 
## Node number 12: 4473 observations,    complexity param=0.01118905
##   mean=65.7087, MSE=5461.028 
##   left son=24 (1898 obs) right son=25 (2575 obs)
##   Primary splits:
##       Iws   < 14.975 to the right, improve=0.11437580, (0 missing)
##       month < 3.5    to the right, improve=0.05671232, (0 missing)
##       DEWP  < 1.5    to the left,  improve=0.03522381, (0 missing)
##       TEMP  < 15.5   to the right, improve=0.02571441, (0 missing)
##       hour  < 7.5    to the right, improve=0.01468811, (0 missing)
##   Surrogate splits:
##       DEWP  < 0.5    to the left,  agree=0.725, adj=0.352, (0 split)
##       hour  < 9.5    to the right, agree=0.636, adj=0.143, (0 split)
##       month < 5.5    to the left,  agree=0.601, adj=0.060, (0 split)
##       TEMP  < 29.5   to the right, agree=0.579, adj=0.008, (0 split)
##       Ir    < 8.5    to the right, agree=0.577, adj=0.003, (0 split)
## 
## Node number 13: 14291 observations,    complexity param=0.02122454
##   mean=110.9471, MSE=6245.967 
##   left son=26 (12080 obs) right son=27 (2211 obs)
##   Primary splits:
##       month < 9.5    to the left,  improve=0.028517730, (0 missing)
##       DEWP  < 24.5   to the left,  improve=0.019623580, (0 missing)
##       TEMP  < 17.5   to the right, improve=0.015208470, (0 missing)
##       Iws   < 3.125  to the right, improve=0.013368580, (0 missing)
##       day   < 17.5   to the left,  improve=0.009606513, (0 missing)
##   Surrogate splits:
##       PRES < 1019.5 to the left,  agree=0.885, adj=0.26, (0 split)
## 
## Node number 14: 697 observations
##   mean=90.76758, MSE=5315.475 
## 
## Node number 15: 3837 observations,    complexity param=0.01763057
##   mean=182.2541, MSE=14770.37 
##   left son=30 (3339 obs) right son=31 (498 obs)
##   Primary splits:
##       month < 1.5    to the right, improve=0.07767789, (0 missing)
##       DEWP  < -7.5   to the left,  improve=0.03723637, (0 missing)
##       PRES  < 1026.5 to the right, improve=0.02702919, (0 missing)
##       hour  < 13.5   to the left,  improve=0.02630411, (0 missing)
##       TEMP  < -2.5   to the right, improve=0.01944054, (0 missing)
##   Surrogate splits:
##       Is   < 14.5   to the left,  agree=0.872, adj=0.016, (0 split)
##       TEMP < -8.5   to the right, agree=0.872, adj=0.014, (0 split)
## 
## Node number 24: 1898 observations
##   mean=36.59852, MSE=2560.741 
## 
## Node number 25: 2575 observations
##   mean=87.16544, MSE=6513.792 
## 
## Node number 26: 12080 observations,    complexity param=0.01018529
##   mean=105.2373, MSE=4881.66 
##   left son=52 (9882 obs) right son=53 (2198 obs)
##   Primary splits:
##       DEWP    < 21.5   to the left,  improve=0.043127390, (0 missing)
##       month   < 3.5    to the right, improve=0.031511180, (0 missing)
##       day     < 17.5   to the left,  improve=0.015688690, (0 missing)
##       cbwd_NE < 0.5    to the right, improve=0.012635700, (0 missing)
##       TEMP    < 8.5    to the right, improve=0.009504751, (0 missing)
## 
## Node number 27: 2211 observations,    complexity param=0.02122454
##   mean=142.1429, MSE=12548.68 
##   left son=54 (1694 obs) right son=55 (517 obs)
##   Primary splits:
##       DEWP < 10.5   to the left,  improve=0.290281800, (0 missing)
##       PRES < 1023.5 to the right, improve=0.109387400, (0 missing)
##       Iws  < 26.54  to the right, improve=0.040596870, (0 missing)
##       day  < 4.5    to the left,  improve=0.013207720, (0 missing)
##       TEMP < 11.5   to the left,  improve=0.008913221, (0 missing)
##   Surrogate splits:
##       PRES < 1012.5 to the right, agree=0.77, adj=0.017, (0 split)
## 
## Node number 30: 3339 observations
##   mean=169.1728, MSE=12214.46 
## 
## Node number 31: 498 observations
##   mean=269.9618, MSE=23067.33 
## 
## Node number 52: 9882 observations
##   mean=98.39425, MSE=4545.826 
## 
## Node number 53: 2198 observations
##   mean=136.0032, MSE=5234.468 
## 
## Node number 54: 1694 observations
##   mean=108.8005, MSE=8101.065 
## 
## Node number 55: 517 observations
##   mean=251.3926, MSE=11543.56

rsq.rpart(reg_tree)

## 
## Regression tree:
## rpart(formula = train_set$pm2.5 ~ ., data = data.frame(train_set), 
##     method = "anova", control = list(minsplit = 5))
## 
## Variables actually used in tree construction:
## [1] cbwd_NW DEWP    Iws     month   TEMP   
## 
## Root node error: 249697509/30279 = 8246.6
## 
## n= 30279 
## 
##          CP nsplit rel error  xerror     xstd
## 1  0.089483      0   1.00000 1.00004 0.015193
## 2  0.067676      1   0.91052 0.91295 0.014715
## 3  0.027921      2   0.84284 0.84944 0.013428
## 4  0.021225      3   0.81492 0.82017 0.013465
## 5  0.019772      5   0.77247 0.76936 0.012906
## 6  0.017631      6   0.75270 0.75775 0.012755
## 7  0.014697      7   0.73507 0.74148 0.012205
## 8  0.011189      8   0.72037 0.73405 0.012073
## 9  0.010185      9   0.70918 0.71550 0.011924
## 10 0.010000     10   0.69900 0.70959 0.011900

Como podemos observar en las gráficas anteriores la regresión por este método no es demasiado buena, ya que con el set de entrenamiento en el mejor de los casos tenemos un valor de R² bastante bajo de alrededor de .3 y errores relativos bastante altos superiores a .7. A pesar de lo anterior, buscaremos probar nuestro modelo con el set de entrenamiento con un numero de splits de 7 y por tanto un CP de .014697. Lo anterior es debido a que vemos que a partir de 7 splits, el error relativo y el R² ya no cambian demasiado.

Ahora veamos una gráfica de los residuales sobre el set de validación:

pruned_tree <- prune(reg_tree, cp=.014697)
pruned_tree_predictions_val <- predict(pruned_tree, data.frame(validation_set))
residuals_pruned_tree_val <- validation_set$pm2.5 - pruned_tree_predictions_val
df_2_val <- data.frame(fitted_values=pruned_tree_predictions_val,
                 residuals=residuals_pruned_tree_val)


ggplot(df_2_val, aes(x=fitted_values, y=residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color="red") +
  theme_bw()

MSE <- sum(residuals_pruned_tree_val^2)
R2 <- 1 - sum((validation_set$pm2.5-pruned_tree_predictions_val)^2)/sum((validation_set$pm2.5-mean(validation_set$pm2.5))^2)

Como podemos observar de la gráfica anterior tenemos residuales relativamente altos, ya que vemos que hay bastantes puntos sobrepasando un residuo de 200, además, hay un MSE de 1.984569110^{7} y un MSRE de 4454.8502955. Lo anterior, junto con que hay un R² de 0.2983416 nos dice que este modelo no es muy bueno para predecir los niveles de PM2.5.

Como hemos visto en esta sección, un árbol de regresión no es lo suficientemente bueno para lo que queremos ya que tiene un R² bastante bajo, sin embargo, en la siguiente sección buscaremos utilizar muchos árboles de regresión para hacer una regresión mucho mejor, es decir, que utilizaremos Random Forest para hacer una regresión.

Random Forest

Debido a que vimos que el rendimiento de un solo árbol de regresión era algo bajo, buscaremos mejorar la regresión con una gran cantidad de árboles. Utilizando distintos parámetros de numero de árboles (100, 200, y 300), hemos visto que 200 dan un muy buen rendimiento, ya que tampoco queremos hacer overfitting del modelo. A su vez, intentamos eliminar la variable de dirección del aire, es decir, la variable cbwd, sin embargo, el rendimiento del modelo se vio disminuido.

rf_reg <- randomForest(pm2.5 ~., data=data.frame(train_set), 
                       ntree=200,
                       nodesize=7,
                       importance=T)
varImpPlot(rf_reg)

Con la gráfica anterior podemos ver que las variables menos importancia tienen para este modelo son las variables de Ir, IS y Cbwd. Por el contrario la variable de temperatura de rocío, día, mes y velocidad del aire son de las más importantes, así como la presión y la temperatura. Como notábamos en el Data Profiling y con el análisis exploratorio vemos que el día y el mes tienen una importancia en el modelo.

A continuación veamos la gráfica de los residuales:

fitted_values_rf1_val <- predict(rf_reg, 
                         newdata=data.frame(validation_set),
                         type="response")

residuals_rf1_val <- validation_set$pm2.5 - fitted_values_rf1_val


df_rf1_performance_val <- data.frame(fitted_values=fitted_values_rf1_val,
                                residuals=residuals_rf1_val)

df_rf1_performance_val[which(is.na(df_rf1_performance_val$fitted_values)), 
                  "fitted_values"] <- 0

ggplot(df_rf1_performance_val, aes(x=fitted_values, y=residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color="red") +
  theme_bw()

MSE <- mean(residuals_rf1_val^2)
MSE

## [1] 2852.86

MRSE <- sqrt(MSE)
MRSE

## [1] 53.41217

y <- validation_set$pm2.5
R2 <- 1 - sum((y-fitted_values_rf1_val)^2)/sum((y-mean(y))^2)
R2

## [1] 0.7162671

Como podemos ver en la gráfica de los residuales, la mayoría de estos no exceden las 200 unidades, y al mismo tiempo, están mayormente concentrados en residuos menores a 100 unidades, lo cual es una gran mejora respecto a la regresión lineal y a la regresión por un solo árbol. Así mismo, podemos observar que este es el modelo que mejor rendimiento tiene ya que tenemos una R² próxima a 1 de 0.7162671, así como un MSE de 2852.8603127 y un MSRE de 53.4121738. Debido a los motivos anteriores este será el modelo que utilizaremos para nuestro set de pruebas.

Mejor modelo en el set de pruebas

Hagamos nuestra predicción sobre el set de pruebas:

fitted_values_rf1_val <- predict(rf_reg, 
                         newdata=data.frame(test_set),
                         type="response")

residuals_rf1_val <- test_set$pm2.5 - fitted_values_rf1_val


df_rf1_performance_val <- data.frame(fitted_values=fitted_values_rf1_val,
                                residuals=residuals_rf1_val)

df_rf1_performance_val[which(is.na(df_rf1_performance_val$fitted_values)), 
                  "fitted_values"] <- 0

ggplot(df_rf1_performance_val, aes(x=fitted_values, y=residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color="red") +
  theme_bw()

MSE <- mean(residuals_rf1_val^2)
MSE

## [1] 4171.298

MRSE <- sqrt(MSE)
MRSE

## [1] 64.58559

y <- test_set$pm2.5
R2 <- 1 - sum((y-fitted_values_rf1_val)^2)/sum((y-mean(y))^2)
R2

## [1] 0.5231877

Podemos observar que la mayoría de los residuales se encuentran alrededor del 0, lo cual nos da indicios de que nuestro modelo tiene un rendimiento bastante decente. Además, los residuales se encuentran en su mayoría por debajo de las 100 unidades lo cual es bastante razonable tomando en cuenta la desviación estándar de todo el set de datos es de aproximadamente 92 unidades. Así mismo tenemos un R² de 0.5231877, así como un MSE de 4171.2981801 y un MSRE de 64.585588.

Conclusiones

A lo largo de este trabajo pudimos observar que las variables descritas en el set de datos para nada tenían un comportamiento lineal respecto a la variable que indicaba los niveles de contaminación de Beijing. Es por ello que se recurrió a utilizar árboles para mejorar las predicciones, y definitivamente pudimos mejorarlas.

Examen_licenciatura_137926_148792

Álvaro Casillas y Ricardo Lara

12/6/2017